﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "lrgs")]
    public static unsafe double lrgs(double a, double b, double eps, IntPtr f_x_ptr)
    {
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return lrgs(a, b, eps);
    }

    /// <summary>
    /// 勒让德_高斯求积法
    /// f计算被积函数值f(x)的函数名。
    /// </summary>
    /// <param name="a">积分下限。</param>
    /// <param name="b">积分上限。要求b>a。</param>
    /// <param name="eps">积分精度要求</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double lrgs(double a, double b, double eps = 1.0E-6)
    {
        int m, i, j;
        double s, p, ep, h, aa, bb, w, x, g = 0.0;
        double[] t = new double[5] {
            -0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459
        };
        double[] c = new double[5] {
            0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851
        };

        m = 1;
        h = b - a;
        s = Math.Abs(0.001 * h);
        p = 1.0e+35;
        ep = eps + 1.0;
        while ((ep >= eps) && (Math.Abs(h) > s))
        {
            g = 0.0;
            for (i = 1; i <= m; i++)
            {
                aa = a + (i - 1.0) * h;
                bb = a + i * h;
                w = 0.0;
                for (j = 0; j <= 4; j++)
                {
                    x = ((bb - aa) * t[j] + (bb + aa)) / 2.0;
                    w = w + f_x(x) * c[j];
                }
                g = g + w;
            }
            g = g * h / 2.0;
            ep = Math.Abs(g - p) / (1.0 + Math.Abs(g));
            p = g;
            m = m + 1;
            h = (b - a) / m;
        }
        return (g);
    }
    /*
    // 勒让德_高斯求积法例
      int main()
      { 
          double a,b,eps,g,lrgsf(double);
          a=2.5; b=8.4; eps=0.000001;
          g=lrgs(a,b,eps,lrgsf);
          cout <<"g = " <<g <<endl;
          return 0;
      }
    // 计算被积函数值
      double lrgsf(double x)
      { 
          double y;
          y=x*x+sin(x);
          return(y);
      }
    */
}

